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ABSTRACT 

\ 

It  is  known  chat  the  stable  GI/FH/1  queue  has  an  eabedded  Markov  chain 
whose  invariant  probability  vector  is  matrix-geometric  with  a  rate 
matrix  R.  In  terms  of  the  matrix  R,  the  stationary  waiting  time  distribu¬ 
tions  at  arrivals,  at  an  arbitrary  time  point  and  until  the  customer's 
departure  may  be  evaluated  by  solving  finite,  highly  structured  systems  of 
linear  differential  equations  with  constant  coefficients.  Asymptotic  results, 
useful  in  truncating  the  computations,  are  also  obtained.  The  queue  is 
assumed  to  follow  the  first-come,  first-served  discipline. 

/ 
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L.  The  Algorithmic  Procedure 

This  paper  is  a  sequel  Co  [2].  All  notations  and  definitions,  which 
were  introduced  there,  will  also  be  used  here  and  have  the  sane  meaning. 

The  purpose  of  the  present  paper  is  to  show  that,  for  the  GI/PH/1  queue, 
the  stationary  probability  distributions  of  the  waiting  time  of  an  arriving 
customer  and  of  the  virtual  waiting  time,  as  well  as  that  of  the  tlme-in- 
system  may  be  computed  by  the  numerical  integration  of  highly  structured, 
finite  systems  of  differential  equations  with  constant  coefficients. 

Only  waiting  times  under  the  first-come,  first-served  discipline  are  considered 
These  systems  of  differential  equations  require  prior  evaluation  of  the 
matrix  R,  which  is  shown  in  [2]  to  be  the  minimal,  nonnegative  solution  to 
a  nonlinear  matrix  equation  of  the  form 

(1)  R  -  I  RkAt  . 

k«*Q  k 

In  the  stable  GI/PH/1  queue,  the  matrix  R  is  positive  and  its  Perron  eigenvalue 
n  satisfies  0<  n  <  1.  The  probabilistic  significance  of  the  matrix  R  is 
discussed  In  [3]  and  a  detailed,  general  treatment  of  matrix-geometric 
solutions  in  stochastic  models  may  be  found  in  [5]. 


A.  The  Waiting  Time  at  Arrivals 


Is  given  by 


(2)  W*(s)  -  C  o  l  Rn  [(sI-T)"  T°A°  J1 e,  for  Re  s  >  0  . 

n-0  ” 

From  Formula  (2) ,  the  mean  waiting  time  W  may  be  computed  by  routine 
differentiations.  One  obtains 

(3)  W  -Xj*  l2  -  , 

where  L2  is  the  mean  queue  length  at  an  arbitrary  time.  The  quantity  L£  is 
explicitly  given  as  a  function  of  R  and  the  parameters  of  the  queue  by 
Formula  (92)  of  [2].  Upon  rewriting  (3)  as  L2  -  X^-1  (WfW^' ) ,  we  recognize 
the  classical  relation,  known  as  Little^  Formula. 

He  now  consider  the  matrix 
00 

(4)  **(s)  -  l  Rn  f (sI-T) -1T0A0 Jn  , 

n-0 

whose  entries  are  Laplace-Stieltjes  transforms  of  mass-functions  on  [0,«) .  The 
matrix  * *(s)  clearly  satisfies 

(5)  **(a)  -  I  +  R  **(s)(8l-T)-1T0A0  . 

In  order  to  obtain  the  equation  (5)  in  a  more  transparent  form,  ve  consider 
the  vector  $*(s) ,  obtained  by  forming  the  direct  sum  of  the  rows  of  4*(s).  The 
vector  v  is  similarly  obtained  from  the  identity  matrix  I. 

Equation  (5)  may  then  be  equivalently  written  as 

(6)  4>*(s)  -  v  +  ^*(s)[R^*(sI-T)-Xt0A0], 


where  <g) denotes  the  Kronecker  product  of  matrices  and  RT  is  the  transpose  of  R. 

The  latter  equation  may  now  be  transformed,  by  using  classical  properties 
of  the  Kronecker  product,  into  the  following  successive  forms: 

&*(b)  -  v+i*(s)  [i®  <sI-T)_1]  (Rt®T°A°)  , 

4*(s)  [I®  (sI-T)-1]^®  I-I®)T  -  RT(g)T°A0]  -  v  , 
and  finally 

(7)  £*(s)  -  v  +  v[sl  ®I-I®T-Rt®  TOA0]-1^^  T°A°) 

The  existence  of  the  matrix  Inverse  in  (7)  will  be  proved  below. 

Let  <K«)  be  the  m2-vector  of  mass-functions,  corresponding  to  the  vector 
of  Laplace-Stleltjes  transforms  J*(s).  It  then  readily  follows  from  (7)  that 

(8)  <Kx)  -  v  +  v  f  exp  [(I®T  +  RT®T°A®>u]du(RT®T°A0) 

“  J0 

-  v  +  v  (I<2)T  +  RT®T°A0)-1  {exp  [(I®T  +  RT<X>T°A°)x] 

-  I(xU}  (RTg)T°A°)  , 

for  x  >  0. 

If  we  Introduce  the  m2-vectors  v°  and  jKx)  by  setting 

(9)  y°  -  -v(I(gT  +  Rt®T°A°)_1  , 

and 

(10)  e(x)  -  v°  exp  [<I®T  +  RT®T°A°)x]  ,  for  x>  0 


then  Formula  (8)  leads  to 


(11)  £(x)  -  v  +  v°(RT(?)T0A0)  -  £(x)  (RT(2)  T°A°)  ,  for  x  >  0. 

Let  now  ♦  (•)»  9(0  end  V°  be  the  mxm  matrices  of  which  the  direct  sums  of 

2  *  a 

the  rows  are  respectively  the  ■  -'vectors  £(  •) »  j9(»)  and  v  .  Formulas 

(11)  and  (9)  may  then  be  rewritten  as 

(12)  *(x)  -  I  +  R  V°T°A0-  R  9  (x)  T°A°  ,  for  x  >  0  , 

and 

(13)  V°T  +  RV°T°A°  -  -  1  . 

From  Equation  (10),  it  is  clear  that  the  vector  £(*)  satisfies  the 
matrix-differential  equation 

JB '  (x)  -  8(x)(I<S>T  +  RT<S>T°A°),  £(0)  -  v°  , 

or  equivalently 

(14)  e'(x)  -  9(x)T  +R0(x)  T°A°  ,  0(0)  -  V°. 

The  matrix  V°  may  be  explicitly  expressed  in  terms  of  R.  To  see  this, 
we  postmultlply  in  (13)  by  the  column  vector  e_  and  recall  that  Tj;  -  -T°.  He 
then  obtain 

V°T°  -  (I-R)-1e  , 

so  that  V°T°A°  may  be  replaced  by  (I-R)“*e  •  a  ■  (I-R)~^A°°.  Doing  so,  readily 


yields 


MSS  - 


MM 
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(15)  V°  -  (I-R) -1(R-I-RA°°)T-1  . 

Equation  (13)  corresponds  to  an  Inhomogeneous  system  of  m2  linear  equations 
In  m2  unknowns.  The  coefficient  matrix  I&T  +  RT &  T°A°  of  the  latter  is 
now  clearly  nonsingular.  A  more  detailed  result  is  proved  in  Lemma  2  below. 
From  Formulas  (2)  and  (12) ,  we  now  obtain  that 

W(x)  -  C  +  C  a  R  V°T°  -  C  a  R  6(x)T°  ,  for  x  >  0. 

Since  however,  a  RV°T°-  -  a  R(I-R) -1e  ,  and  C  -  [a/I-R)-1e]-1  , 

this  expression  may  be  further  simplified  to 

(16)  W(x)  -  1-  C  a  R  0(x)T°  ,  for  x  >  0  . 

The  results  of  this  derivation  may  be  summarized  as  follows. 

Theorem  1 

If  for  the  stable  GI/PH/1  queue,  the  matrix  R  is  known,  then  the  stationary 
waiting  time  distribution  W(*)  is  given  by  Formula  (16).  The  matrix  0(x)  is 
obtained  by  numerical  integration  of  the  matrix-differential  equation  (14). 

Remark 

We  note  that  Formula  (16)  provides  a  transparent  solution  to  Lindley's 


equation  for  a  wide  subclass  of  GI/G/1  queues.  For  this  subclass,  the  solution 
is  more  elementary  and  appropriate  for  numerical  implementation  than  the 
classical  approach,  which  is  based  on  Wiener-Hopf  techniques. 
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The  stationary  distribution  W^( •)  of  the  time-in-system  of  an  arriving 
customer  may  be  easily  computed  along  with  the  probability  distribution 
H(*).  One  appropriate  computational  organization  is  discussed  in  Section  4.1 
of  [5];  we  present  a  somewhat  different  one  here. 

The  distribution  W^(*)  is  the  convolution  of  W(*)  and  the  service  time 
distribution  H(*>.  Since  the  latter  is  of  phase  type  with  Irreducible 
representation  (<*,T> ,  we  have 

(17)  W  (x)  *  W(u)f*  exp[T(x-u)]T°du  ,  for  x  >  0  . 

1  J0 

If  we  set 

(18)  £(x)  “of  tf(u)exp  [T(x-u)]  du  ,  for  x  ^  0  , 

J  o 

then  elementary  calculations  show  that 

(19)  £'(x)  -  £(x)T  +  W(x)o  ,  for  x  >  0  , 

with  £(0)~  0..  In  order  to  compute  W^(*)»  it  therefore,  suffices  to  integrate 
the  system  of  differential  equations  (19)  and  to  evaluate  W^(x)  “  £(x)T°,  for 
x  >  0.  We  note  that  this  is  a  general  procedure  for  the  numerical  convolution 
of  two  distributions  on  [0,“) ,  when  one  of  them  is  of  phase  type. 

B.  The  Virtual  Waiting  Time 

In  Theorem  8  of  [2],  the  Laplace-Stleltjes  transform  W*(s)  of  the  virtual 

A 

waiting  time  distribution  W( •)  is  given  by 

m 

(20)  W*(s)  -  1  -  p  +  l  CoRi~1  f[R][(sI-T)"1T0A°]ie,  for  ReB  >_  0  . 

i-1 


f 


The  matrix  T[R]  is  defined  in  Formula  (74)  of  [2]  and  satisfies  the 
equation 

(21)  R  -  I  +  X[  1>[R]  +  Xx'  R?[R]  T°A°  , 

proved  in  Lemma  9  of  the  same  paper.  The  matrix  f  [R]  may  be  explicitly  obtained 
in  terms  of  R.  This  was  overlooked  in  [2]  and  the  following  lemma  Indeed 
yields  the  result  of  Theorem  6  more  directly. 

Lemma  1 

The  matrix  f [r]  is  given  by 

(22)  T[R]  -  X '  “1(R-I-R  A00)T-1. 

Proof 

Postmultiplication  by  e  in  (21)  yields 

Re  -  e  -  X^  f[R]T°  +  Xj/R  $  [R]  T°. 

Since  I-R  is  nonsingular,  it  follows  that  X^'  T[r]  T°  “  £,  and  hence 

Xi'  y [R]  T°A°  •  A°°.  Upon  substitution  into  (21),  we  readily  obtain  Formula  (22). 

The  transformation  of  Formula  (20) ,  required  to  obtain  a  convenient 
algorithm  for  the  computation  of  the  probability  distribution  H(*),  is  entirely 
analogous  to  that  given  above.  We  shall  only  show  the  most  important  steps. 

We  may  write  (20)  as 

(23)  W*(s)  -  1-p  +  C  b  S*(s)(sI-T)"1T0  , 

where  •*(»)  is  anm  x  m  matrix  of  Laplace~Stieltjes  transforms,  satisfying 


(24) 


♦*(s)  -  (R-I-RA00) T-1  +  R  $*(8)(8l-T)"1T°A° 


Performing  the  same  manipulations  as  those  following  Equation  (5) ,  we 
obtain 

(25)  fcx)  “  (R-I-RA°°)T-1  +  RV°T°A°  -  R  0  (x)T°A°  , 
for  x  ^  0.  This  equation  corresponds  to  (12)  above. 

A 

The  matrix  V°,  which  now  satisfies  the  equation 

(26)  V°T  +  R  V°T°A°  -  -(R-I-RA00)T-1  , 
is  explicitly  given  by 

(27)  V°  -  - (R-I-RA00) T-2  -  R(I-R)"1(R-I-RA00)T-1A00T“1 

-  -v°T-l +  RV°(I-A00)T"1  , 
where  V°  is  given  in  Formula  (15)  . 

The  matrix  0 (•)  is  obtained  by  solving  the  matrix-differential  equation 

(28)  0*(x)  ■  9(x)T  +  R  0  (x)T°A°  ,  0(0)  -  V°. 

A 

As  we  compute  the  matrix  ♦  (•)»  by  use  of  (28)  and  (25),  we  simultaneously 
integrate  the  differential  equation 

(29)  3 (x)  -  ^'(x)!  +  i (x),  fj<0)  -  0. 

This  yields  the  inverse  of  the  Laplace-Stieltjes  transform  **(s) (sI-T)-^-. 

A 

The  probability  distribution  W( •)  is  then  finally  given  by 


The  algorithm  for  W ( * )  is  now  clear.  We  evaluate  the  matrix  V°  and  inte¬ 
grate  the  differential  equations  (28)  and  (29)  by  classical  methods.  These 
will  usually  give  the  computed  value  of  4^(x)  at  values  of  x,  which  are  the 
successive  multiples  of  an  appropriately  chosen  step  h.  We  note  that  it  is 

A  A 

not  necessary  to  store  all  the  preceding  values  of  *(•)  and  #^(»),  but 

only  those  few  that  are  needed  by  the  integration  procedure.  For  each  computed 

A  A 

value  of  #^(*)»  the  corresponding  value  of  W(*)  is  evaluated  and  stored 
for  printout. 

A 

The  mean  of  W(«)  is  obtained  by  a  routine  differentiation  in  Formula  (20) 
and  is  given  by 

(31)  W  -  p  W  +  \  \i'2 

where  vi 2'  is  the  second  moment  of  the  service  time.  The  quantity  W  is  the 
stationary  mean  waiting  time  at  arrivals. 

2  .  Asymptotic  Results  for  the  Waiting  Time  Distributions 

The  preceding  derivations  also  readily  yield  precise  asymptotic  results 

A 

for  the  probability  distributions  W(»),  W^(’)  and  W(*).  We  first  discuss 
some  preliminary  matters. 

The  positive  matrix  R  has  the  maximal  eigenvalue  q,  satisfying  0  <  q  <  1. 
Let  the  vectors  and  u°  be  left  and  right  eigenvectors,  corresponding  to  q. 
Both  vectors  may  be  chosen  to  be  positive  and  to  satisfy  u  e  ■  u  u°  «■  1. 

The  matrix  T  +  q  T°A°  is  clearly  an  irreducible,  stable  matrix.  It 


therefore  has  an  eigenvalue  -£  <  0,  which  is  simple  and  for  which  the 
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corresponding  left  and  right  eigenvectors  z  and  z°  may  be  chosen  to  be  positive 
and  to  satisfy  z_  ±  m  z.  z°  ■  1.  Any  other  eigenvalue  £'  satisfies  Re(£')  <  -E, 
The  eigenvalue  -£  will  be  called  the  eigenvalue  with  maximum  real  part  of 
T  +  n  T°A°. 


Lemma  2 

a.  The  vectors  u  and  z_  are  identical. 

b.  The  matrix  I(S)T  +  RT ®  T°A°  has  nonnegative  off-diagonal  elements. 

Its  eigenvalues  £"  satisfy  Re(5*')  <_  -£,  and  the  quantity  -£  is  its 
eigenvalue  of  maximum  real  part.  The  left  and  right  eigenvectors 
corresponding  to-£  are  given  by  u°  (S)u_  and  u T(g>£°,  respectively.  The 
inner  product  of  these  two  positive  vectors  is  equal  to  one. 

Proof 

We  recall  from  [2]  that,  in  the  case  of  the  GI/PH/1  queue  ,  the  matrix 

oo 

A*(z)  -  £  A  zv,  0  <  z  <  1,  is  given  by 

v*0  v 

(32)  A*(z)  -  (  exp[(T+zT°A°)t]  dF(t) . 

•'0 

Moreover,  q  is  also  the  maximal  eigenvalue  of  the  positive  matrix  A*(n)  and 
the  vector  u  is  the  corresponding  left  eigenvector,  whose  components  sum 
to  one. 

On  the  other  hand,  it  readily  follows  from  (32)  that  z_  A*(n)  "  f*(£)z.* 
where  f*(‘)  denotes  the  Laplace-Stleltjes  transform  of  F(«).  The  uniqueness 
of  the  Perron  eigenvalue  and  of  the  (normalized)  corresponding  left 
eigenvectors  now  readily  imply  that  u  ■  £,  and  q  ■  f*(£). 


■  ..  • 

&  V 
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It  Is  clear,  from  Che  definition  of  Che  matrix  Ig)T  +  R‘  <§)T°A°,  that 
Its  off-diagonal  elements  are  nonnegative.  Moreover,  the  positivity  of  R 
and  the  irreducibility  of  T+T°A°  Imply  that  this  matrix  is  also  irreducible. 

We  have 

(u°T®  u)  (I  ®  T  +  RT  ®  T°A°)  -  u°T®  uT  +  u^R1®  uT°A°  - 
u°<guT  +  r\u°® uT°A°  -  u°T<g)u(T  +  nT°A°)  -  -£(u°(g>u)  . 

X 

A  similar  calculation  holds  for  the  vector  u^x)  z°.  Since  the  vectors  u°®u 
and  a*  positive,  they  are  respectively  the  left  and  right  eigenvectors 

corresponding  to  the  eigenvalue  -£  of  maximum  real  part  of  the  matrix 
I  (x)  T  +  Rt  (g>  T°A° .  Finally  (u0*^  u>  (ulg)  £°)  -  (u  u°)T®  (u  z}  “  !• 

It  now  follows  from  Formula  (10)  that 

(33)  0(x)  -  -v(I®T  +  Rt (x) T°A°) ~X [(u^  z°)  •  (u°fe>  u)  ]e~^x  +  o  (e"tx)  , 
a  8  x  -*■  °».  The  coefficient  of  e~^*  may  be  simplified  as  follows 

-v (I(x)  T  +  RT®> T°A°)  _1  [(wfe £°)  •  (u°^5>  u)] 

-  i  v  [(u%)20).(u°^u)] 

£ 

-_i.  (u  z°)(u°^)u)  -  i  (u°^g)u)  . 

£  "  “ 

T 

The  m^-vector  u°®u  is  the  direct  sum  of  the  rows  of  the  mxm  matrix  u.°u,  so 
that  Formula  (33)  may  be  equivalently  written  as 

(34)  0(x)  -  i  u°  u  e“^x  +  o(e“£x)  ,  as  x  +  •. 

6 
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Theorem  2 

The  probability  distribution  W(*)  satisfies 

(35)  l-W(x)  -  k  e“^x  +  o(e“^x),  as  x  +  •. 

The  constant  k  is  given  by 

(36)  k  -  Cn(l-n)_1(o  U°). 

Proof 

Upon  substitution  of  (34)  into  (16) ,  we  obtain 

l-W(x)  *  £  a  R(u°u)T°  e“^x  +  o(e“^x)  ,  as  x  -*■  ••  . 

However,  Ru°  "  nu°,  and  uT  +  n(uT°)a  “  -£u.  Upon  postmultlpllcatlon  by_e 
in  the  latter  equation,  we  obtain  that  u  T°  ■  C(l-n)-1.  The  atated  result 
is  now  immediate. 

It  is  also  clear  from  the  last  two  equations  in  the  preceding  proof 
that  u  -  -£n (1-n) a (£I+T) .  Upon  postmultlpllcatlon  by  T°,  we  obtain  the 
useful  formula 

(37)  -aUt+T)-1  T°  -  I  . 

n 

Since  W^(*)  satisfies  the  equation  (17),  one  readily  obtains  that 

(38)  1-W-^x)  -  -a(CI+T)T°  •  ke"*x  +  0(e“*x) 

-  n'^k  e~^x  +  o(e"^x),  as  k  -*•  •  . 


•t'  >■  ■  * 


An  entirely  similar  argument  holds  for  the  distribution  W(*).  One 

obtains  after  some  calculations  that 

(39)  l-W(x)  -  k  e"5x  +  o(e“5x)  , 

vhere  k  is  given  by 

(40)  k  -  -A^_1  £_1(l-n)a(Cl+T)_1I°  •  k 

-  A^"1  r1  C(o  u°)  . 

Remarks 

a .  The  asymptotic  results  in  Formulas  (35) ,  (38)  and  (39)  correspond  to 
the  exact  results  for  the  Gl/M/1  queue.  For  that  simple  case,  the 
term  o(e~^x)  may  be  omitted;  the  resulting  formulas  then  hold  for 
all  x  ^  0. 

b.  He  shall  show  below  how  the  quantities  q  and  £  may  be  computed  by  an 
elementary  and  efficient  algorithm.  We  then  see  from  the  proof  of 
Theorem  2  that 

(41)  u  -  -Cn(l-n)"1  a(£l+T)_1  . 

This  formula  provides  us  with  a  powerful  accuracy  check  in  the  computation 
of  the  matrix  R.  The  latter  is  obtained  by  iterative  solution  of 
Equation  (1) .  The  quality  of  the  computed  value  for  R  may  be  assessed 
by  seeing  how  closely  the  equation  uR  “  qu,  is  satisfied.  The  amount  of 
computation  required  to  obtain  £,q  and  ti  is  small,  compared  to  that 
needed  to  get  R. 


e.  It  is  Impossible,  in  gensral,  to  compute  the  vector  u°  without 

evaluating  the  matrix  R.  The  coefficients  k,  k  and  k,  therefore, 
require  a  fair  amount  of  computation.  Even  without  precise  knowledge 
of  these  coefficients,  the  asymptotic  results  have  a  number  of  practical 
uses.  We  propose  to  discuss  these  in  a  forthcoming  paper,  which  deals 
with  similar  asymptotic  results  of  much  greater  generality. 

The  following  result  provides  us  with  an  elementary  numerical  procedure, 
which  simultaneously  yields  the  values  of  n  and  £.  We  note  that  the  Laplace- 
Stleltjes  transform  of  the  service  time  distribution  is  given  by 
h(a)  “  a(sI-T)“*T°.  This  transform  is  a  rational  function  in  s  and  has  an 
abscissa  of  convergence  — r*  <  0.  The  function  h(s)  is  strictly  decreasing 
on  the  interval  (-• r*,«0  and  satisfiea  h(0)  -  1. 

We  define  ifr(z)  to  be  the  unique  real  solution  of  the  equation 

(42)  o[p(z)I-T]_1T0  -  i  , 

where  0  <  z  <_  1.  We  see  that  p(z)  increases  from  -t*  to  zero,  as  z  varies 
from  0  to  1. 

Theorem  3 

For  the  stable  GI/PH/1  queue,  the  quantity  n  is  the  unique  solution  in 
(0,1)  to  the  equations 

(43)  *  -  f*[-*(*)j  ,  I-hO(z)J  , 
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Proof 

Tho  quantity  f (z)  la  tha  eigenvalue  of  max lmm  real  part  of  the  marrlx 
T+sT°A°.  To  see  this,  we  write 

u(z)  [T+  z  T°A°]  -  p(z)u(z)  ,  0  <  z  *  1  , 

and  normalize  u(z)  by  u(z)jK  “  1. 

Ue  then  readily  obtain 

u(z)T°  -  -Ij.  p(z)  , 

and  also 

u<*)  »  *(z)o  [♦(z)I-t]"1  . 

Pos  multiplication  by  T°  yields  (42) .  It  Is  also  elementary  to  verify  that 
the  vector  u(z)  Is  positive.  The  quantity  f(z)  Is  therefore  the  eigenvalue 
of  maximum  real  part  of  the  stable  matrix  T+zT°A°,  for  0  <  s  £  1. 

It  Is  now  lamadiate  that 

u(z)A*(z)  -  u(z)  [  exp[(T+zT°A°)u]  dF(u)  -  f*[-q>(z)]  u(z)  , 

so  that  f*[-<i(z)]  Is  the  Perron  eigenvalue  of  the  positive  matrix  A*(s)  for 
0  <  z  <  1. 

Ue  know  however  from  [2]  that,  provided  the  queue  is  stable,  n  Is 
the  unique  solution  in  (0,1)  of  the  equation  z  ■  f*[-f(*)]>  Finally,  from 
the  definition  of  ,  It  is  clear  that  £  ■  H*(n). 


J 
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To  compute  n  and  £,  we  nay  either  solve  the  equation  f*(-s)h(e)  “  1, 
for  Its  unique  solution  s°  in  (-t*,0)  and  set  £  ■  -sq ,  and  n  ”  f*(C)  “  or 
we  may  compute  i|i(z)  for  successive  values  of  z  by  using  (42)  and  select  the 
latter  to  approximate  the  solution  to  z  ■  f*[-f(z)]  in  (0,1).  Any  of  the 
classical  methods,  such  as  the  bisection,  secant,  or  Newton's  method  may 
be  Implemented.  Convergence  is  rapid  and  the  required  computation  time  is 
very  small. 

3.  Asymptotic  Behavior  of  the  Quene  Length  Densities 

The  powers  of  the  matrix  R  have  the  well-known  asymptotic  behavior 

( 45)  ■  nk  u°  •  u  +  o  (n^) ,  as  k  -*•  •  . 

Furthermore,  it  follows  from  (22)  that 

u  *[R]  -  1^_1  [(n-Du  -  no  ]  T"1  . 

Using  the  expression  for  u  obtained  in  (41),  we  may  routinely  simplify  the 
preceding  formula  to 

(46)  u  iMR]  -  -Xj"1  na(5l4T)"1. 

The  following  results  are  now  immediate. 

Theorem  4 


In  the  stable  GI/PH/1  queue,  the  stationary  queue  length  densities  satisfy 


19 


(47)  -  C  a  Rk  -  C(a  u°)  nku  +  o(nk) 
jgfcg.  “  C(a  u°)qk  +  o(nk), 

l  2Ejj  *  C(a  u°) (1-n)-1  ?jk  +  o(nk)  » 

v-k 

(48)  -  Ca  Rk_1,f[R]  -  -X|  _1  C(a  u°)  nk  •  o(Cl-KT)”1  +  o(nk)  , 
lit  S."  _xi  1  c<2.  “°>  nk  o(Cl+T)_1e  +  o(nk)  , 

ee 

X  Zk  “  -xi_1  c^sl  H°)  *  nk  •  i(c1+r)-x  +  o(nk) 

v-k 

(49)  *k  “  [  I  y  I0]"1  Ifc+i  T°  -  CoRke  -  c(a  u°)  nk  +  o(nk)  , 

v-1  ^ 

m 

l  \  -  c(a  u0)(l-n)_1  nk  -h>(nk)  , 

V-k 

as  k  ■*  m. 

Formulae  47-49  give  asymptotic  expressions  for  the  stationary  densities 
and  dlatrlbutlons  of  the  queue  length,  respectively  at  arrival  epochs,  at  an 
arbitrary  time  and  immediately  after  departure  epochs. 

Remarks  and  Acknowledgments 

In  Section  5  of  [2],  ve  described  a  recursive  procedure  to  compute  the 
matrices  A^,  k  0,  for  the  PH/PH/1  quene.  This  procedure  Is  particularly 
appealing  as  It  completely  avoids  the  use  of  numerical  Integration.  Ve  did 
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express  concern,  however,  over  the  possible  singularity  of  the  matrix 
1(5) S  +  I(gl,  which  appears  in  it.  This  concern  is  without  foundation. 

As  the  Kronecker  sum  of  the  nonsingular  matrices  S  and  T,  the  matrix  I(g)S  + 
T®I  is  nonsingular  [i] .  The  recursive  procedure  for  the  computation  of  the 
matrices  Ak  has  now  been  Implemented  and  performs  very  well  [4] . 

The  author  is  indebted  to  Professors  Y.  Takahashl  and  V.  Ramaswaml 
for  helpful  comments  on  the  derivations  in  this  paper. 
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